home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / histogram / stat.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  3.2 KB  |  141 lines
  1. /* gsl_histogram_stat.c
  2.  * Copyright (C) 2000  Simone Piccardi
  3.  *
  4.  * This library is free software; you can redistribute it and/or
  5.  * modify it under the terms of the GNU General Public License as
  6.  * published by the Free Software Foundation; either version 2 of the
  7.  * License, or (at your option) any later version.
  8.  *
  9.  * This program is distributed in the hope that it will be useful,
  10.  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12.  * General Public License for more details.
  13.  *
  14.  * You should have received a copy of the GNU General Public
  15.  * License along with this library; if not, write to the
  16.  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  17.  * Boston, MA 02111-1307, USA.
  18.  */
  19. /***************************************************************
  20.  *
  21.  * File gsl_histogram_stat.c: 
  22.  * Routines for statisticalcomputations on histograms. 
  23.  * Need GSL library and header.
  24.  * Contains the routines:
  25.  * gsl_histogram_mean    compute histogram mean
  26.  * gsl_histogram_sigma   compute histogram sigma
  27.  *
  28.  * Author: S. Piccardi
  29.  * Jan. 2000
  30.  *
  31.  ***************************************************************/
  32. #include <config.h>
  33. #include <stdlib.h>
  34. #include <math.h>
  35. #include <gsl/gsl_errno.h>
  36. #include <gsl/gsl_histogram.h>
  37.  
  38. /* FIXME: We skip negative values in the histogram h->bin[i] < 0,
  39.    since those correspond to negative weights (BJG) */
  40.  
  41. double
  42. gsl_histogram_mean (const gsl_histogram * h)
  43. {
  44.   const size_t n = h->n;
  45.   size_t i;
  46.  
  47.   /* Compute the bin-weighted arithmetic mean M of a histogram using the
  48.      recurrence relation
  49.  
  50.      M(n) = M(n-1) + (x[n] - M(n-1)) (w(n)/(W(n-1) + w(n))) 
  51.      W(n) = W(n-1) + w(n)
  52.  
  53.    */
  54.  
  55.   long double wmean = 0;
  56.   long double W = 0;
  57.  
  58.   for (i = 0; i < n; i++)
  59.     {
  60.       double xi = (h->range[i + 1] + h->range[i]) / 2;
  61.       double wi = h->bin[i];
  62.  
  63.       if (wi > 0)
  64.         {
  65.           W += wi;
  66.           wmean += (xi - wmean) * (wi / W);
  67.         }
  68.     }
  69.  
  70.   return wmean;
  71. }
  72.  
  73. double
  74. gsl_histogram_sigma (const gsl_histogram * h)
  75. {
  76.   const size_t n = h->n;
  77.   size_t i;
  78.  
  79.   long double wvariance = 0 ;
  80.   long double wmean = 0;
  81.   long double W = 0;
  82.  
  83.   /* FIXME: should use a single pass formula here, as given in
  84.      N.J.Higham 'Accuracy and Stability of Numerical Methods', p.12 */
  85.  
  86.   /* Compute the mean */
  87.  
  88.   for (i = 0; i < n; i++)
  89.     {
  90.       double xi = (h->range[i + 1] + h->range[i]) / 2;
  91.       double wi = h->bin[i];
  92.  
  93.       if (wi > 0)
  94.         {
  95.           W += wi;
  96.           wmean += (xi - wmean) * (wi / W);
  97.         }
  98.     }
  99.  
  100.   /* Compute the variance */
  101.  
  102.   W = 0.0;
  103.  
  104.   for (i = 0; i < n; i++)
  105.     {
  106.       double xi = ((h->range[i + 1]) + (h->range[i])) / 2;
  107.       double wi = h->bin[i];
  108.  
  109.       if (wi > 0) {
  110.         const long double delta = (xi - wmean);
  111.         W += wi ;
  112.         wvariance += (delta * delta - wvariance) * (wi / W);
  113.       }
  114.     }
  115.  
  116.   {
  117.     double sigma = sqrt (wvariance) ;
  118.     return sigma;
  119.   }
  120. }
  121.  
  122.  
  123. /*
  124.   sum up all bins of histogram
  125.  */
  126.  
  127. double 
  128. gsl_histogram_sum(const gsl_histogram * h)
  129. {
  130.   double sum=0;
  131.   size_t i=0;
  132.   size_t n;
  133.   n=h->n;
  134.  
  135.   while(i < n)
  136.     sum += h->bin[i++];
  137.  
  138.   return sum;
  139. }
  140.  
  141.